Localization in one dimensional lattices with non- nearest- neighbor hopping: 
Generalized Anderson and Aubry- Andre models 
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We study the quantum localization phenomena of noninteracting particles in one-dimensional 
lattices based on tight-binding models with various forms of hopping terms beyond the nearest 
neighbor, which are generalizations of the famous Aubry- Andre and noninteracting Anderson model. 
For the case with deterministic disordered potential induced by a secondary incommensurate lattice 
(i.e. the Aubry- Andre model), we identify a class of self-dual models, for which the boundary 
between localized and extended eigenstates are determined analytically by employing a generalized 
Aubry- Andre transformation. We also numerically investigate the localization properties of nondual 
models with next-nearest-neighbor hopping, Gaussian, and power-law decay hopping terms. We 
find that even for these non-dual models, the numerically obtained mobility edges can be well 
approximated by the analytically obtained condition for localization transition in the self dual 
models, as long as the decay of the hopping rate with respect to distance is sufficiently fast. For 
the disordered potential with genuinely random character, we examine scenarios with next-nearest- 
neighbor hopping, exponential, Gaussian, and power-law decay hopping terms numerically. We find 
that the higher-order hopping terms can remove the symmetry in the localization length about the 
energy band center compared to the Anderson model. Furthermore, our results demonstrate that 
for the power-law decay case, there exists a critical exponent below which mobility edges can be 
found. Our theoretical results could, in principle, be directly tested in shallow atomic optical lattice 
systems enabling non- nearest-neighbor hopping. 
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I. INTRODUCTION 

Quantum transport of matter waves in the presence of 
disorder has long been a topic of interest for condensed- 
matter physicists. For one-dimensional (ID) noninter- 
acting systems, one of the oldest and most extensively 
studied models for quantum transport is the single-band, 
nearest-neighbor (nn) tight-binding model, 



ti{u n +i + u n -i) + V n u n = Eu n , 



(1) 



where t\ is the nn hopping integral term and V n is the 
on-site disordered potential 1 . One of the main merits of 
Eq. ([!) is its simple form, which lends itself to fast nu- 
merical analysis, as well as exact theoretical statements 
on quantum transport in some cases. Arguably the most 
well known of the latter is when V n is bounded, uncorre- 
cted disorder (i.e., the noninteracting Anderson modeP), 
where it can be shown that all eigenstates of the system 
are localized for any nonzero potential strength. Another 
well-known example is the ID incommensurate problem 
[in particular, V n = V cos(27ran + J), where a is irra- 
tional] studied by Aubry and Andre where all eigenstates 
of the system are extended for potential strength below 
a threshold value (V c = 2ti) and localized above this 
threshold^. Conversely, this simple form given by Eq. 
([!]) makes it difficult to study directly in solid-state sys- 
tems, where the disorder is difficult to control reliably 
and interactions can rarely be ignored. However, recent 
advances in the manipulation of ultracold atoms in op- 
tical lattices provide a powerful tool for directly exam- 



ining quantum transport in fundamental models such as 
Eq. 0. This has been most notably demonstrated in re- 
cent experiments conducted by Billy et al., where Ander- 
son localization was directly observed in a Bose-Einstein 
condensate (BEC) subjected to a random laser speckle 
potentiaP, and similarly in experiments conducted by 
Roati et a/., where Aubry- Andre duality was directly 
observed in a BEC loaded into a quasi-periodic optical 
lattice^. These feats, which previously eluded experi- 
mental observation for decades, illustrate the potential 
of ultracold atomic systems to experimentally probe fun- 
damental quantum phenomena. Considering the degree 
of control afforded in ultracold atomic systems, we can 
systematically relax, in a controlled manner, basic as- 
sumptions inherent in Eq. ([!]) and directly study their 
influence on quantum transport and how it differs from 
the well-known Anderson and Aubry- Andre results. It is 
this potential in ultracold atomic systems that motivates 
our present work, where we examine quantum transport 
in tight-binding, non-interacting models that are exten- 
sions of Eq. 0. In particular, we relax the nn tight- 
binding assumption and theoretically examine transport 
in models with long and short-range hopping schemes. 
Such models should be representative of diffuse gases of 
ultracold atoms loaded into fairly shallow optical lattices. 



We can go beyond the nn coupling assumption while 
still remaining in the tight-binding framework by includ- 
ing higher-order hopping terms. The general form of such 
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a model with on-site disorder is given by 



(2) 



where the tight-biding terms t m may assume a variety 
of forms. There is a small, but growing body of nu- 
merical and analytical work examining transport in the 
context of this generalized modeP^. In an effort to 
extend this growing body of work, we wish to investi- 
gate transport through various forms of Eq. Q with 
both incommensurate and random on-site potentials. In 
particular, we study transport in tight-binding models 
with next-nearest-neighbor (nnn) hopping (the t\ — 
model) and models incorporating an infinite number of 
hopping terms that decay by an exponential, Gaussian, 
or inverse power law. Since we examine both deter- 
ministic bichromatic potentials and random potentials, 
this report is divided into two sections. The first sec- 
tion examines the case of incommensurate potentials by 
first studying the exponential hopping model, which has 
been shown to have an analytical mobility edge^, then 
approximately extrapolating these results to predict the 
mobility edges in the nnn, Gaussian, and inverse-power- 
law hopping models^. The second section numerically 
examines Eq. ([2| with nnn, exponential, Gaussian, and 
inverse-power-law hopping for randomly disordered po- 
tentials and highlights how localization in these models 
is markedly different compared to what is seen in the case 
of the nn Anderson model. 



II. INCOMMENSURATE POTENTIALS 

A. Self Dual models 

One of the first models to examine quantum transport 
in ID incommensurate potentials is the so-called Aubry- 
Andre (AA) model 3 . In this model, the on-site term 
m Eq. is a cosine with frequency incommensurate 
with the primary lattice: V n = V cos (2ir an + S) where 
a is an irrational number and S is an arbitrary phase. 
It has been shown that this model is self-dual under the 
transformation, 



im(27ran+5) ifin 



(3) 



when V = 2t\ (i.e. f m satisfies the same eigenvalue equa- 
tion as u n ). Also, under Eq. ([3]), if the eigenstate u n is 
spatially localized, then the eigenstate of the dual prob- 
lem, / m , is spatially extended and vice versa. Using this 
property and the Thouless formula for incommensurate 
potentials 16 , it is argued that all eigenstates are localized 
for V > 2ti and extended for V < 2t\. The case where 
V = 2t\ is especially interesting and has been shown to 
yield a singular continuous eigenspectrum, which forms 
a Cantor set in the thermodynamic limit 17 . Furthermore 



AA duality can be shown to have a more general forirPR 
Consider the model, 



y^JmUn+m + V V n U n Eu n . 



(4) 



If the on-site potential and the hopping terms satisfy the 
relation, 

V n =Y, t mZ im{2 ™ n+S \ ( 5 ) 



then this model also possesses an AA-like duality. 

Other models have be en sh own to possess self-duality 
similar to the AA modeJMLyj. The particular model we 
consider here, given by 

Eu n = te- p \ n - n '\u n > + Ucos(27ran + S)u n , (6) 

is especially interesting because its self-duality condition 
naturally predicts energy-dependent mobility edges, in 
contrast to the AA modef^H To see this we now reiterate 
the results given in the previous work on this model to 
show that Eq. (|6| does have a self-dual condition and we 
expand on the previous work to show that this self-dual 
condition does indeed define a mobility edge. 

We begin by defining the parameter, po > such that 

(E + t)-V cos(27ran + S) = ftT n , (7) 
cosh(po) - cos(27ran + S) 

1 n — r-j-7 \ , (p) 

smh(po) 

n = V(E + t) 2 -V 2 . (9) 

Then it follows that (E + t)/V = cosh(po) and we can 
rewrite Eq. (|6| as 



Y^ te -p\n-n'\. 



If we now consider the transformation 

~ _ im{2Kan+5)r r 
771 — / j - L n Uj n- l 

n 

and note that for p > we have the identity 



rn 1 = E e " p 



\m\ im(27ran+5) 



(10) 



(11) 



(12) 



then it follows that the state u m satisfies the equation 
nf m u m = J2 te- po|m - m,| fw, (13) 



where T m is given by 

^ cosh(p) — cos(27rara + S) 
7,1 sinh(p) 



(14) 



We see that Eq. (10) is self-dual under the transforma- 



tion Eq. (11) when p = po, or equivalently cosh(p) 



3 



cosh(po) for p,p > 0- Thus, the duality condition for 
Eq. (|6| is given by 



cosh(p) 



E + t 
V ' 



(15) 



Since the transformation given by Eq. (11) transforms 



localized states to extended states and vice versa [sim- 
ilar to Eq. ([3|], we expect that the eigenstates of the 
system are critical (neither localized nor extended) when 
Eq. (15) is satisfied. 



Similar to the arguments made by Aubry and Andre 
for the AA model, we now expand on the conjecture made 
in the previous workP^I and argue that the eigenstates of 
Eq. ( 10 ) are localized for p > po and extended for p < po 
[i.e. that Eq. (15) does, indeed, define a mobility edge]. 



Since the Thouless formula used by Aubry and Andre 
was derived for models of the form of Eq. ([!]), we can 
not use it for our particular model. Therefore, our first 
step is to generalize the idea of the Thouless formula for 
the long-range hopping model. To do so, we treat Q as 
the eigenvalue and consider the Green's matrix, 



G(Q) m , n = (01 -H) 
cofactor (£11 — H) 

ILCn-fy) 



1 



(16) 
(17) 



where the cofactor is the appropriately signed determi- 
nant with the mth row and nth column removed and 
H is the Hamiltonian corresponding to the eigenvalue 
equation given in Eq. ( [To]) where we have set t = 1 
without loss of generality; I is the identity matrix. As- 
suming a nondegenerate eigenspectrum, the Green's ma- 
trix has a simple pole for each eigenvalue, Q^. Since, 
by definition, the residue of G(0^) m?n is the product 
of the mth and nth elements of the eigenvector (i.e. 

Res(G(fi,) m ,n, ^/i) = ^m^n^jP^, we have for the product 
of the first and last elements of the eigenvector 



cofactor^^I — H) 



1,7V 



rw^ - sip) 



(18) 



If the state is exponentially localized about the site n', 
then we expect u n ~ exp(— j\n f — n|), where 7 > is 
the characteristic (or Lyapunov) exponent. Therefore the 
product uiUn ~ exp(—j(N — 1)). Thus, the characteris- 
tic exponent for large N is given by 



7 (fi M ) = lim -(N 

N^oo 



l)- 1 ln|u^ ) «^ ) l 



N 



lim (N — l) -1 (y^ In \Slfj, — Slp\ 

- In I cofactor (^1 - H)i jJV |). (19) 

This is the generalized Thouless relation for the charac- 
teristic exponent of a wavefunction. For the case where 
H is given by Eq. (10), the cofactor takes on the form: 



cofactor(^ /i I — H) 



1,7V 



N-2„-(N-l)prp-l 

1 N ' 



(20) 



Then we have for the characteristic exponent, 
7(n M )=p-ln|n^|+ lim (jV-l)-iy ln|O p 

/V — Vr«n ' ^ 



Sir- 



(21) 



We now compare the characteristic exponents of the 
eigenvectors of Eq. (10), which we denote as 7(0), with 
the exponents of the dual problem Eq. (13), denoted as 



7(0). Since the eigenvalue Q is not changed by the trans- 
formation given by Eq. (11), we expect the summation 



term on the right-hand side of Eq. (21) to be equal for 



both Eq. (10)andEq. (13). Therefore, the characteristic 



exponents have the following relation: 

7(H)=7(fi) + (p-Po). 



(22) 



Considering the case when p > po, since > 0, it 

follows that 7(^) > and therefore the eigenstate, u n is 
localized while the dual state, u n is extended. Similarly, 
when po > p, we can argue that 7(f2) > and there- 
fore the dual state, u n , is localized while u n is extended. 
Therefore, returning to the original problem given by Eq. 
([6| and using the fact that cosh(po) is a monotonically 
increasing function of po > 0, it follows that the eigen- 
states are localized for (E + t)/V < cosh(p) and extended 
for (E + i)/V > cosh(p). 

Similar to the AA model, the self-duality described 
above has a general form, considering, again, a model of 
the form given in Eq. Q . The model will have this form 
of self-duality if the on-site potential and the hopping 
terms satisfy the relation: 



B 



E 



j. im(27ran+5) 



(23) 



where A and B are constants. In particular, the constant 
B gives the slope of the the duality condition (i.e., B = 
E/V ). 

We now numerically examine localization in Eq. ([6| 
[and equivalent ly Eq. (10)] by calculating the inverse 
participation ratio (IPR) of the wave functions, given by: 



IPR 1 



EnK 



W|4 



(24) 



where the superscript i denote the i eigenstate. The 
IPR approaches zero for spatially extended wave func- 
tions and is finite for localized wave functions, and hence 
has a useful diagnostic role. Fig. [l] plots energy eigen- 
values (or eigenstate number) and the IPR of the cor- 
responding wave functions for Eq. (|6| as a function of 
potential strength, V, with a = — l)/2 and p = 1, 
2, 3, or 4. The solid curves in the figures represent the 
boundary given in Eq. (15). From the figure we see that 



IPR values are approximately zero for energies above the 
boundary and are finite for energies below the bound- 
ary. This supports our assertion that the mobility edge 
is, indeed, given by Eq. (15). 
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In Fig. [2j we directly examine sample eigenstates in 
each regime (i.e., localized, extended, and near the mo- 
bility edge) for p = 1.5 and V = 1.8. We see that the 
wavefunction is localized for low energies [Fig [2^l] , ex- 
tended for high energies [Fig [2Jd] , and critical near the 
boundary [Fig [2^]. 



Finally, we examine the eigenvalues of Eq. ( 10) for dif- 
ferent values of a at the duality point (p = popwheie we 
expect the eigenspectrum to form a fractal set for large 
TV. The results of this are given in Fig. [3j In the fi sure, 
we see that for large values of the eigenspectrum closely 
resembles the well-known Hofstadter's butterfly whi ch re- 
sults from the solutions of Harper's equatio n 119 * 20 *. For 
smaller values of however, we see a generalized form of 
Hofstadter's butterfly that is not symmetrical about the 



the band center, but skewed toward lower eigenvalues. 
The self-similarity in the figure suggests that the eigen- 
spectrum does, indeed, form a Cantor set at the duality 
point in the thermodynamic limit. 



B. Nondual models 

General realizations of Eq. Q with an incommensu- 
rate potential should not be expected to satisfy either 
Eq. ([5| or Eq. (23). Thus, the mobility edges in non- 
dual incommensurate problems may not be discernible 
by theoretical means. However, approximate theoretical 
statements can be made for some non-dual models with 
hopping terms that fall off in some general manner. If 
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FIG. 2. Eigenstates of Eq. ([6} with 500 lattice sites, 
a = (y/E — l)/2, V = 1.8, and p = 1.5 for different energy 
eigenvalues: a) low energy localized state below the mobility 
edge b) high energy extended state above the mobility edge 
c) critical state near the mobility edge. 



this fall off is fast enough, then the localization transi- 
tions in these models are largely determined by the ratio 
W^i- This ratio can be used to determine an approxi- 
mately equivalent model of the form of Eq. |6]), which, 
as shown in the above section, has an exact theoreti- 
cal localization boundary. To show this possibility, we 
numerically examine the eigenstates of several non-dual 
models. In particular, we examine the t\ — model, 
the Gaussian hopping model, and the inverse-power-law 
hopping model with an incommensurate potential and 
see how closely the approximate localization boundary 
matches with the numerically observed one. 

The t\ —t<i model is the nnn extension of the AA model 
and is given by 



+ U n -2) + h(u n +i + U n -i) 

+ V cos(27ran + S)u n 



Eu n . (25) 



The parameters of the approximately equivalent expo- 
nential hopping model [Eq. (|6|] are given by p = 
ln(ti/t2) and t = t\e p '. Using Eq. (15), we can ap- 



proximate the boundary between localized and extended 
states. To examine this approximation we calculate the 
IPR of the eigenstates of Eq. ( |25| ) . The results are given 
in Fig. [4] for 500 lattice sites, a = (\/5 — l)/2 and vari- 
ous values of the ratio £2/^1 ■ The solid lines in the figure 
give the approximate mobility edge given by Eq. (15). 



From the figure, we see that for small values of t^Jti, 
the approximate boundary is in good qualitative agree- 
ment with the numerical IPR results. For larger values 
(h/h > 0.3), however the boundary differs considerably 
from the linear condition in Eq. (15). 




FIG. 3. Eigenspectr um o f Eq. (10) with varying a for (a) 
p = 1 , (b) p = 2, and (c) p = 3. 



The tight-binding incommensurate model with Gaus- 
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FIG. 4. Inverse participation ratios and energy eigenvalues of all eigenstates of Eq. ( 25 ) with 500 lattice sites and a — (V§ — 1)/2 



for £2/^1 — |(a)| 0.05, (b) 0.1, (c) 0.3, |(d)| 0.5. The solid curves represent the approximate boundary between spatially localized 
and spatially extended states. 



sian hopping has the form: 



^2 e -°\ n - n ' 



V cos(27ran + 5)u n = Eu n . (26) 



Similar to the t\ —£2 model, the approximately equivalent 
exponential hopping model can be determined from the 
ratio £2/^1, which yields p = 3<r. The IPR results for 
this model are given in Fig. [5] [again, a = (s/E — l)/2]. 
In this figure, we see that the approximate boundary is 
in good qualitative agreement with the numerical results 
for larger values of a. Small values of a, however, result 
in very interesting energy-dependent mobility edges that 
are not linear in potential strength, which is similar to 



the £1 — £2 results for large values of £i/£2- 

For power law decay in the hopping terms, we examine 
the model: 



V cos(27ran + S) = Eu n 



(27) 



In this case, the exponential coefficient, p is given by 
p = rln(2). Fig [6] gives IPR results for this model with 

2 [Fig 



a = (Vb - I)/ 2 and r 1/3 [Fig 6(a 



|6(b)| , and r = 3 [Fig |6(c)] . In each of these cases, the 
approximate localization boundary is in good qualitative 
agreement with the numerical results. 

From the above numerical results, we see that the lo- 
calization boundary for the exponential hopping model 
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FIG. 5. Inverse participation ratios and energy eigenvalues of all eigenstates of Eq. (26) with 4096 lattice sites and a — 
0/5- l)/2 for a = (a) 1.0, (b) 0.5, (c) 0.25, and (d) 0.10. The solid curves represent the approximate boundary between 
spatially localized and spatially extended states. 



gives a good qualitative agreement for the Gaussian and 
inverse-power-law hopping models with large-enough de- 
cay coefficients (i.e., r and a > 1). Thus, we believe that 
the localization boundary for any tight-binding model 
with hopping terms that decay fast enough can be ap- 
proximated by the results of the exponential hopping 
model. In general, however, we see that the energy- 
dependent mobility edges in the non-dual models are not 
linear in V as with the case in the exponential hopping 
model. Whether there is an exact theoretical statement 
to describe these peculiar mobility edges is still an open 
question. 



III. RANDOM DISORDER 

We examine the case of a random potential in the tight 
binding framework in the context of both the nnn model 
(i.e. the t\ — 1<2 model) and extended schemes where the 
coupling may be short ranged in the sense of decaying 
exponentially (or more rapidly, as in a Gaussian decay) 
or long-ranged with a power law decay. 

We examine the effects of a random potential directly 
in the context of characteristics of the eigenstates, by 
studying the IPR, which provides information as to the 
extent of localization. We calculate the IPR versus eigen- 
state number for a range of random potential strengths, 
which will tend to localize states. 

We produce surface graphs of the IPR to show the 
characteristics of the eigenstates with respect to localiza- 
tion. In addition, one may calculate histograms of the 
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(c)r = 3 



FIG. 6. Inverse participation ratios and energy eigenvalues of all eigenstates of Eq. ( |27[ ) with 2000 lattice sites and a = 
(y/E — l)/2 for various values of r. The solid curves represent the approximate boundary between spatially localized and 
spatially extended states. 



IPR instead of preparing surface plots with respect to 
eigenstate number and the width of the potential. The 
former complement the latter by showing how the sta- 
tistical weight for a particular IPR value evolves with 
increasing system size. 

We find in both the surface plots and in the histograms 
a convergence toward bulk behavior, where self- averaging 
in sufficiently large systems reduces the differences in the 
characteristics of the eigenstates corresponding to inde- 
pendently generated random potential realizations. In 
the thermodynamic limit, the type of tight-binding cou- 
pling scheme and the overall statistical characteristics of 
the random potential are the factors which determine the 
distribution of the properties of the eigenstates. 

In implementing the on-site random potential, we op- 



erate either in terms of a binary potential where the po- 
tential at a particular site assumes the value V or —V 
with equal probability or a continuous one which is cho- 
sen with uniform probability between the bounds V and 
— V. In either case, the strength of the random potential 
may be considered to be parameterized by V. We begin 
by studying the t\ — t<i model with random on-site dis- 
order. Fig. [7] gives the IPR as a function of the disorder 
strength, V, and the eigenstate number, z, for the eigen- 
states of the t\ — t<i tight-binding model with a binary 
random on-site potential and 5000 lattice sites. In the 
figure we see that for =0, we have the ID Anderson 
model where localization is relatively weak at the band 
center and comparatively strong near the band edges 1 . 
However as we increase t<i, we see that the weakly local- 
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ized states shift to higher energies toward the top band 
edge. Moreover, as the relative value of ti is increased, 
cross sections of the IPR along the eigenstate number 
axis develo p a bimodal profile. In the case where £ 2 = 0.8 
[Fig. |7(d)| the bimodal structure is particularly stark. In 
this case we begin to see two distinct regions of weakly 
localized states, in contrast to the Anderson case where 
the localization is weakest in the band center, giving way 
to more strongly localized states at the band edges. We 
also obtain qualitatively similar results when we consider 
a uniformly distributed, rather than binary, random po- 
tential (Fig. ||). 

In addition to surface plots produced for a single real- 
ization of the random potential, we also generate surface 
plots of the IPR with respect to V and eigenstate number 
where the results are averaged over a large number of ran- 
dom configurations. In particular, in Fig. [9j we average 



the IPR over 100 different realizations of an on-site, uni- 
formly distributed random potential for the t\ — ti model 
with 500 lattice sites. The qualitative similarities among 
Figures [7] [8] and [9] suggest that the localization behavior 
is self-averaging in the sense that statistical fluctuations 
play a small role in determining the characteristics of the 
system if the system size TV is sufficiently large. 

We now study the randomly disordered tight-binding 
model with a Gaussian decay in the hopping terms. 



Fig. [10] and Fig. [TT] give the IPR as a function of po- 
tential strength, V, and eigenstate number, z, of the 
eigenstates of the Gaussian decay tight-binding model 
[i.e. tm = exp(— a|m| 2 )] with a random binary and ran- 
dom uniform on-site potential, respectively, and 5000 lat- 
tice sites. In both Fig. [To] and Fig. [TT] we see that for 
small a, the weakly localized states appear at higher en- 
ergies in comparison to the Anderson model (similar to 



10 




1.5 2 2.5 3 3.5 4 




1.5 2 2.5 3 3.5 4 



(a)t 2 = 



(a)t 2 = 




1 1.5 2 2.5 3 3.5 4 

V 

(b)t 2 = 0.2 

FIG. 8. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the t\ — £2 model with a random, uniformly 
distributed potential for 5000 lattice sites. 



the t\ — £2 model for large t^). For larger a, the IPR 
is qualitatively similar to that of the Anderson model 
case as would be expected given the very rapid decay 
which strongly suppresses hopping beyond the nearest 
neighbors for which the Anderson model is an idealiza- 
tion with hopping confined strictly to nearest neighbors. 
We also average the IPR over 100 different realizations 
of an on-site, uniformly distributed random potential for 
the Gaussian decay hopping model with 500 lattice sites 
and report the results in Fig. [12) A gain, the qualitative 
similarities among Fig. [TUJ Fig. |11[ and Fig. [12] suggest 
that self- averaging is at work in the localization char- 
acteristics of the eigenstates. The averaging over many 
realizations of disorder has the effect of removing much of 
the graininess due to statistical fluctuations which would 
not appear in the bulk limit and seem to be finite size 
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(b)t 2 = 0.2 

FIG. 9. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the t\ — £2 model averaged over 100 realizations 
of random, uniformly distributed potentials with 500 lattice 
sites. 



artifacts, while preserving more smoothly varying char- 
acteristics which appear to be associated with the bulk 
limit. 

Similar results can be seen in the randomly disordered 
tight-binding model with exponentially decaying hopping 
terms. Figures 13 and [14] give the IPR as a function 
of potential strength, V, and eigenstate number, z, of 
the eigenstates of the exponential hopping tight-binding 
model [i.e. t m = exp(— p|ra|)] with a random, binary 
and random, uniform on-site potential, respectively, and 
5000 lattice sites. Similar to the Gaussian decay model, 
we see that for small p, the weakly localized states appear 
at higher energies compared to the Anderson case, and 
for larger p, the IPR approaches that in the Anderson 
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FIG. 10. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the Gaussian hopping model with a random, 
binary potential for 5000 lattice sites. 



FIG. 11. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the Gaussian hopping model with a random, 
uniformly distributed potential for 5000 lattice sites. 



case. Moreover, just as we saw in the t\ — and Gaus- 
sian models, evidence of self- averaging may be seen by 
examining the average of the IPR over 100 realizations 
of a uniform random potential with 500 lattice sites (Fig. 
15). 

We now turn to the randomly disordered tight-binding 
extended model with hopping terms that decay by an 
inverse- power- law. Since this slow form of decay essen- 
tially allows for long-range hopping, the localization char- 
acteristics of this model differ more from those of the An- 
derson model, where the hopping scheme is short-ranged, 
than those of the other models we have investigated. In 



Figures 16 and[TTJ we show the IPR for a random, binary 
and random, uniform on-site potential respectively as a 
function of potential strength, V, and eigenstate num- 
ber, z, of the eigenstates of the tight-binding model with 



hopping terms that fall off as t m = l/\m\ r . In these fig- 
ures, we see that the states near the top band edge are 
weakly localized, and the top band edge is possibly de- 
localized for r = 1/3 and r = 1. This supports earlier 
theoretical and numerical work which suggests that there 
is, indeed, a mobility edge at the top band edge for long- 
range hopping^MI]. We revisit this question by examining 
the statistical distribution of the IPR for successive sys- 
tem size doublings. As in the case of the previous models 
with short-range hopping, we again find evidence of self- 
averaging in the surface IPR plots in this model with hop- 
ping terms decaying slowly, as power laws, by examining 
the average of the IPR over 100 realizations of a uni- 
form random potential with 500 lattice sites (Fig. 18). 



Again, averaging over multiple realizations of disorder 
yields a smoother IPR plot by removing minor noisy fea- 
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FIG. 12. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the Gaussian hopping model averaged over 100 
realizations of random, uniformly distributed potentials with 
500 lattice sites. 



tures which are essentially statistical fluctuations about 
the bulk limit. 

The surface plots show the structure of the inverse par- 
ticipation ratio with respect to measures of the random 
potential strength such as V and the eigenstate num- 
ber index. Alternatively, one may concentrate on the 
frequency of particular values of the IPR as a way to 
obtain a statistical description of the characteristics of 
the eigenstates with respect to localization. Although 
we lose specific information for individual eigenstates, we 
gain in return the ability to observe trends in the statis- 
tical characteristics of the IPR distributions with respect 
to increasing system size; in this way, we determine in a 
rigorous manner what portion of the states are localized 
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FIG. 13. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the exponential hopping model with a random, 
binary potential for 5000 lattice sites. 



and what portion, if any, have extended character. 

Localized states are associated with a finite IPR value, 
whereas the IPR will tend to zero for extended states. 
Hence, if in the bulk limit (i.e., in the limit of very 
large N) all of the states for a particular random po- 
tential strength V are localized, the histogram will cease 
to change and take the form of a constant profile indepen- 
dent of TV, which is determined only by the strength of 
the random disorder and the extended coupling scheme. 

On the other hand, if a finite portion of the states have 
genuine extended character, the IPR will continue to de- 
crease for a fraction of the eigenstates, and a portion 
of the histogram total weight will move steadily toward 
lower IPR values. In the case of short ranged models 
such as the nn Anderson model, we find the IPR his- 
togram to eventually shift to a profile which is converged 
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FIG. 14. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the exponential hopping model with a random, 
uniformly distributed potential for 5000 lattice sites. 



with respect to increasing system size TV, and no further 
variation is seen in the shape or position of the IPR his- 
togram curve. With this in mind, we focus our efforts 
on the long-range hopping model (i.e., power-law decay 
in hopping terms) and examine the scaling of the IPR 
histograms with system size to determine the presence of 
mobility edges, which are suggested by the IPR surface 
plots. 

In preparing the IPR histograms, it is important to 
average away statistical fluctuations by sampling a suf- 
ficient number of eigenvalues; for each system size we 
considered, we sample at least n e i g = 4 x 10 5 eigenval- 
ues. We obtain the required statistics by diagonalizing 



samp 



r jN Hamiltonian matrices; the number of 
separate matrices to be considered decreases with in- 
creasing system size N as self-averaging within an in- 
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FIG. 15. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the exponential hopping model averaged over 
100 realizations of random, uniformly distributed potentials 
with 500 lattice sites. 



dividual random potential configuration supplies more 
statistics for larger systems. Hence, whereas 40, 000 ma- 
trices are analyzed for N = 100, for the largest system, 
N = 6400, we examine only 80 random potential realiza- 
tions. 

In Fig. [l9j the tunneling matrix element decreases rel- 
atively rapidly with r = 2.0 for the decay power, and 
the histograms shown in panels (a), (b), (c), and (d) are 
consistent with the scenario in which all eigenstates are 
localized in the bulk limit, even for small values of V or 
weak random potentials. In each of the panels (a), (b), 
(c), and (d) of Fig. 19, histogram curves corresponding to 
various systems sizes N are shown, where the sequence 
of system sizes is chosen to facilitate the study of the 



14 




V 

(a)r = 1/3 




0.5 1 1.5 2 2.5 3 3.5 4 

V 



(b)r = 1 




V 

(c)r = 2 

FIG. 16. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the inverse-power-law model with a random, 
binary potential for 5000 lattice sites. 
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FIG. 17. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the inverse- power- law model with a random, 
uniformly distributed potential for 5000 lattice sites. 
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FIG. 18. (Color online) Surface plots of the base 10 logarithm 
of the IPRs for the inverse-power-law model averaged over 100 
realizations of random, uniformly distributed potentials with 
500 lattice sites. 



effect of successive doubling of N on the histogram pro- 
file. In panels (a) and (b), convergence to an invariant 
histogram curve corresponding to localized states is rel- 
atively swift, while the approach to the limiting profile 
is more gradual in the curves shown in panel (c), where 
V = 0.5, and panel (d) with V = 0.125. Nevertheless, 
the graphs in Fig. [19] indicate a stabilization of the results 
with respect to doubling N for each random potential 
strength shown, and we conclude for the hopping term 
decay exponent r = 2.0 that essentially any random on- 
site potential strength (irrespective of the strength V) is 
sufficient to localize states in the bulk limit). 



For the graphs shown in Fig. 20 , the vertical axis repre- 
sents the base 10 logarithm of the histogram amplitude, 
and the random potential strengths V are identical to 
those of the corresponding graphs in Fig. 19 A salient 



feature of the curves is the convergence to a profile which 
terminates for a particular IPR with no histogram weight 
above this upper limit IPR value. For sufficiently large 
values of the decay exponent r (i.e. for at least r > 2.0), 
there is a minimum value IPR m in where the histogram 
amplitude abruptly falls to zero, and there is no prob- 
ability of finding states with a lower IPR, a condition 
indicating the absolute localization of all states in the 
bulk limit. 

As the decay exponent r is decreased, there seems to 
be a threshold value r c where the properties of the IPR 
histogram change in a qualitative manner with respect to 



increasing N. Fig. 21 and Fig. 22 show histogram curves 
for r = 1.2, whereas Fig. [23] and Fig. [24] display IPR 
histograms for the case r = 1.5 where the decay of the 
tunneling coefficients is faster. A salient common char- 
acteristic in the graphs obtained for r = 1.2 and r = 1.5 
suggests both decay exponents are below the threshold 
value r c . In contrast to the r = 2.0 case, where the his- 
togram profiles converge to a curve which does not change 
with successive doubling of the size TV (a behavior com- 
patible with the localization of all states), for both the 
cases r = 1.5 and r = 1.2, there is a steady advance of 
the leftmost edge of the histogram curve toward smaller 
IPR values in the low-IPR regime. The size of the incre- 
ment appears to be essentially the same each time the 
size of the system is doubled. On the other hand, for 
larger IPR values, the rightmost parts of the histogram 
converge and cease to change with increasing N. 

The shift of the leftward edge toward lower log 10 (IPR) 
values occurs at a constant rate with doubling of the 
size N of the system, a phenomenon seen for all of the 
histograms obtained for r = 1.2 and r = 1.5. While a 
bimodal structure may be seen in both the graphs ob- 
tained for r = 1.2 and r = 1.5, dual peaked character 
is most prominently manifest for the slower decay expo- 
nent r = 1.2 and for lower values of V corresponding to 
weaker random potentials. The two peaks behave very 
differently with increasing N. Whereas the rightmost 
peak, corresponding to relatively higher IPR values and 
hence more localized character does not shift significantly 
in location, the peak on the left migrates steadily toward 
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lower IPR values with successive doubling of N. In addi- 
tion, the peak height appears to decrease at a geometric 
rate each time the system size is doubled. 



The curves shown in Fig. 22 and Fig. 24 , where the 
amplitude of the IPR histogram is presented as a base 
10 logarithm, highlight an important feature for the his- 
tograms in the cases r = 1.2 and r = 1.5 absent in the 
case of the more rapidly decaying scheme where, for ex- 
ample, r = 2.0. For large (but finite) N the log- log 
curves for r < r c may be divided, crudely speaking, 
into three regimes. Moving leftward along the horizon- 
tal log 10 (IPR) axis, one first sees a rapid rise to a maxi- 
mum, and the curve then begins to decrease with decreas- 
ing log 10 (IPR). For intermediate values of the logarithm 
of the Inverse Participation Ratio, the logarithm of the 
histogram amplitude decreases linearly with decreasing 
log 10 (IPR). Finally, the curve rises again to reach a sec- 
ond maximum before beginning a rapid decline. 

The intermediate regime where the logarithm of the 
histogram amplitude decreases linearly is a salient com- 
mon feature which becomes broader as N is increased 
(extending further and further leftward). For r = 1.2 
and V = 0.125, the region where the dependence is ap- 
proximately linear is more difficult to discern, but would 
be more readily seen for systems sizes beyond the largest 
(AT max = 6400) we consider in the context of this study. 
However, even though the linear dependence may not be 
readily visible, the leftward peaks in the log- log graph 
decrease in height at a linear rate even for weaker po- 
tentials (e.g., V = 0.125) where the intermediate linear 
region is more difficult to discern. (The linear decline of 
the peak height in the log-log plot with successive system 
size doublings is compatible with the geometric decline 
evident in the semilogarithmic graphs.) Simple extrapo- 
lation would suggest that as the bulk limit is approached, 
the leftward edge of the histogram curve will continue to 
advance at a constant rate to the left; ultimately, in the 
thermodynamic limit, the asymptotically linear decrease 
in the logarithm of the histogram amplitude would con- 
tinue for arbitrarily small log 10 (IPR). Hence, in terms 
of the histogram density (/>ipr, the support for states de- 
creases with decreasing IPR, ultimately vanishing as the 
IPR heads to zero (i.e., for bona fide extended states). Al- 
though the statistical weighting decreases with decreas- 
ing IPR for all values of r considered, the decline is much 
less abrupt if r < r c . 

In particular, we infer the explicit dependence for the 
IPR density (/>ipr for r < r c will be log lo (0ipR,) = 
a log 10 (IPR) +/?, a relation which would hold in the bulk 
limit for reasonably large values of log 10 (IPR). Inverting 
for IPR yields IPR = /(r, v){IPR) a ^ , where f(r,v) 
and a(r, v) depend on the hopping decay exponent r and 
the random potential strength V (f would also depend 
on the specific type of random disorder, be it generated 
in a binary fashion or sampled from a uniform distribu- 
tion). For fixed potential strength V, we anticipate that 
a(r,v) will rise sharply in the vicinity of r = r c , where 
the decay of the histogram weight becomes much more 
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FIG. 19. (Color online) Histograms are shown versus 
log 10 (IPR) for the decay power r = 2.0 for systems rang- 
ing in size (by successive doubling) from N = 100 to 
N = 6400. Panels (a), (b), (c), and (d) correspond to 
V = 2.0 , 1.0, 0.5 and 0.125, respectively. 



rapid than for r < r c . 

The sudden shift in the behavior of the IPR his- 
tograms as the system size is successively doubled sug- 
gests an abrupt transition from the condition where 
the asymptotic dependence of the histogram density is 
a relatively slow decay with decreasing IPR, 0ipr = 
/(r, V)(IPR) a ^) to a much more rapid decrease. The 
transition is likely signaled by a divergence in the expo- 
nent a at a critical value r c of the decay exponent in the 
extended hopping scheme. 



IV. CONCLUSION 

We have shown with numerical calculations and ana- 
lytical results that (ID) tight-binding models with on-site 
disorder and higher-order hopping terms exhibit interest- 
ing and non-trivial localization phenomena that can vary 
considerably from the well-known results of the nn tight- 
binding models. In the specific case where the on-site po- 
tential is an incommensurate potential (particularly the 
bi-chromatic problem), we have shown that for a gen- 
eral expression for the decay in the hopping terms with 
range, the energy-dependent mobility edges can be pre- 
dicted approximately based on the ratio of the nn and 
nnn hopping terms, txjt^-, for sufficiently fast decay. 

We have also considered the case of bounded, uncor- 
rected disorder, where we have shown that in the exam- 
ined models, the higher order hopping terms can remove 
the symmetry in the localization length about the en- 
ergy band center compared to the nn Anderson model. 
Broadly speaking, it appears that eigenstates with lower 
energies tend to localize with shorter localization lengths 
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FIG. 20. (Color online) Histograms with the amplitude given 
as a base ten logarithm are shown versus log 10 (IPR) for the 
decay power r — 2.0 for systems ranging in size (by successive 
doubling) from N = 100 to N = 6400. Panels (a), (b), (c), 
and (d) correspond to V = 2.0, 1.0, 0.5, and 0.125, respec- 
tively. 



FIG. 22. (Color online) Histograms with the amplitude given 
as a base ten logarithm are shown versus log 10 (IPR) for the 
decay power r = 1.5 for systems varying in size (by successive 
doubling) from N = 100 to N = 6400. Panels (a), (b), (c), 
and (d) correspond to V = 2.0, 1.0, 0.5, and 0.125, respec- 
tively. 
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FIG. 21. (Color online) Histograms are shown versus 
log 10 (IPR) for the decay power r = 1.2 for systems rang- 
ing in size (by successive doubling) from N = 100 to 
N = 6400. Panels (a), (b), (c), and (d) correspond to 
V = 2.0, 1.0, 0.5, and 0.125, respectively. 
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FIG. 23. (Color online) Histograms are plotted versus 
log 10 (IPR) for the decay power r = 1.5 for systems rang- 
ing in size (by successive doubling) from N = 100 to 
N = 6400. Panels (a), (b), (c), and (d) correspond to 
V = 2.0, 1.0, 0.5, and 0.125, respectively. 



(i.e. higher IPR values) compared to higher-energy eigen- 
states due to the presence of the higher-order hopping 
terms. 

There is also the curious case of inverse- power- law hop- 
ping terms where a mobility edge may manifest itself at 
the top band edge if the decay exponent is smaller than 
a certain critical value. 

We have prepared histograms of the IPR to determine 
the statistical characteristics of the IPR distribution for 



the model with power law decay hopping. For relatively 
short-range hopping schemes, the histogram weight falls 
to zero for finite IPR values, suggesting that all states 
are localized in the bulk limit. However, if the hop- 
ping terms decay sufficiently slowly with distance, the 
histogram density is still zero in the low IPR limit, but 
the decay has a much slower asymptotic dependence with 
the form f(r,v)(lPR) a ^ v \ 
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FIG. 24. (Color online) Histograms with the amplitude given 
as a base ten logarithm are shown versus log 10 (IPR) for the 
decay power r = 1.5 for systems ranging in size (by successive 
doubling) from N = 100 to N = 6400. Panels (a), (b), (c), 
and (d) correspond to V = 2.0, 1.0, 0.5, and 0.125, respec- 
tively. 



Our results are especially relevant to current experi- 
mental efforts as we consider alternatives to solid-state 
systems to study quantum transport, such as cold atoms 
in shallow optical lattices. Such systems may not be 
as strongly binding as solid-state systems, and there- 
fore, do not fit well with the nn tight-binding assump- 
tion. Given this and the considerable degree of control 
given to experimentalists in optical lattices, these re- 
sults can be observed experimentally within cold atomic 
systems. In particular, the current experiments in cold 
atomic gases^have already verified the basic (and long- 
established) features of ID quantum localization proper- 
ties in the Anderson and the Aubry model. Experiments 
in shallow lattices allowing longer-range hopping should 
enable a deeper understanding of the localization prop- 
erties discussed in the current work. 

We conclude with some discussion of some of the open 
questions in this subject, which may be of importance for 
future studies. One important issue completely beyond 
the scope of the current work is the effect of interaction 
on our predicted mobility edges in the generalized one- 
particle AA model. In general the solid-state experimen- 
tal systems are many-particle systems, and interaction is 
invariably present. In atomic systems where the inter- 
action is often short-ranged, it is possible to approach 
the noninteracting limit by having a very dilute system, 
and this is one reason behind the recent experimental 
success in studying localization properties in refsF^H It 
is also obvious that, while our one-electron localization 
theory applies both to fermions and bosons since quan- 
tum statistics are irrelevant in the single-particle limit, 
the many-particle interacting situation is different for 
fermions and bosons and must be studied independently. 



In principle this is a formidable task, although in one di- 
mension progress is feasible by combining numerical and 
analytical methods. 

The interacting bosonic problem is easier to study 
theoretically, and for the strict nn- hopping AA model 3 , 
the conclusion — based on extensive density-matrix- 
renormalization-group (DMRG) studies of the problem 
— appears to be that the main features of the minimal 
AA model is preserved although a complex phase dia- 
gram now emerges in the presence of both interaction 
and incommensurate potential manifesting a complex in- 
terplay between superfluid and Mott insulating phases^. 
Such numerical studies using DMRG or perhaps the time- 
evolving-block-decimation (TEBD) technique should, in 
principle, be possible for our extended AA model, and 
it will be an interesting future direction to pursue in 
this problem. An alternative technique for studying the 
bosonic problem in the presence of both the incommen- 
surate AA potential and interaction is to utilize the non- 
linear Schrodinger equation approach using the so-called 
Gross-Pitaevski equation (GPE). Such a study has re- 
cently been carried outP^ for the Anderson model in the 
bosonic case where the interplay of disorder and interac- 
tion was studied both numerically and analytically with 
the conclusion that the basic qualitative feature of the 
non-interacting model is not modified by the presence of 
interaction (i.e. all states remain localized even in the 
presence of interaction). Again, such a GPE-based study 
should, in principle, be possible for our generalized AA 
model, although it is likely to be nu merically challenging. 
Based on the existing results EIE21 our best guess is that 
our conclusion in the current work about the existence of 
a mobility edge in the generalized AA model in the pres- 
ence of non-nearest neighbor hopping will remain valid 
for bosonic systems even in the presence of interaction, 
but more work is necessary to establish this point. 

Studying the interacting system becomes much more 
difficult and complex for fermions where the interplay 
of interaction and localization is notoriously difficult to 
study. It is well known that in general (repulsive) inter- 
action leads to effective derealization since the interact- 
ing particles want to stay away from each other. On the 
other hand, strong interaction also causes Mott transition 
and Wigner crystallization (where the system becomes 
pinned in the presence of disorder), and thus enhances 
localization effects in some situations. Direct numerical 
diagonalization and other works for the fermionic Hub- 
bard model in the presence of disorder^ indicate that 
interaction tends to increase the localization length with- 
out modifying the basic localization properties of the ID 
Anderson model. Although no detailed investigation of 
the AA model for fermions in the presence of interaction 
has yet been undertaken in the literature, it is reason- 
able to assume, based on the results of the correspond- 
ing interacting Anderson model, that the basic picture of 
the noninteracting AA model would remain valid quali- 
tatively even in the presence of interaction. We therefore 
believe that the existence of mobility edges in the general- 
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ized AA model discussed in the current work will remain 
valid even in the presence of interaction for both bosons 
and fermions, but much more work will be needed to set- 
tle this issue decisively. This remains an interesting and 
important problem for future studies. We also mention 
in this context a recent worlP^ which draws an interest- 
ing distinction between the AA model and the Anderson 
model with respect to the nature of the underlying local- 
ization properties, and it will be interesting to investigate 
whether such an analysis sheds insight into our discovery 



of a mobility edge in the generalized AA model in the 
presence of non-nearest-neighbor hopping. 
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